XploR is an R package designed for robust, allelic imbalance and large-scale copy number analysis from whole exome sequencing (WES) data in clinical genomics. It features advanced noise reduction using a panel of normal samples for both coverage and allelic counts, comprehensive smoothing and segmentation algorithms, and accurate purity and ploidy estimation. XploR supports flexible rerun options based on chromosome region, tumor purity, or diploid coverage, and includes integrated ISCN annotation and visualization. These capabilities make XploR a powerful solution for clinical and research applications in genomic copy number analysis.

Contents

• Features
• Installation
• Test Run
• Prepare input files
• Prepare reference files
• Algorithm
• Output
• Model selection and rerun
• Full function and parameter list

Features

• BAF and coverage denoise, smoothing, binning, and quality control
• Exome-wide copy number segmentation and allelic imbalance detection
• Purity and ploidy estimation with model selection
• Rerun based on chromosome region, purity or diploid coverage
• ISCN and gene annotation of CNV segments
• Visualization

Installation

Install the latest version from GitHub using devtools:

install.packages("devtools")
devtools::install_github("sj-cmpb-se/XploR")

Quick Test run

All files needed for a test run in placed at inst/extdata folder. RunExamplePipeline() will use the files in inst/exdata for a test run. Panel of normal generation is not included in the test run. Details for build a panel of normals please refer to Prepare reference files

library(XploR)
RunExamplePipeline( out_dir = "/path_to_output_dir" )

Running this function is same with running the steps separately like:

1. Run segmentation based on Allelic imbalance information. The example used “cbs” segmentation method.

RunAIsegmentation(
    seg = seg,
    cov = cov,
    ai = ai,
    gender = gender,
    out_dir = out_dir,
    prefix = prefix,
    ai_pon = ai_pon,
    aitype = "dragen"
  )

Parameters for RunAIsegmentation

Note on aitype column requirements: - If "gatk" or "other": input must include columns CONTIG, POSITION, ALT_COUNT, REF_COUNT, REF_NUCLEOTIDE, and ALT_NUCLEOTIDE. - If "dragen": input must include columns contig, start, stop, refAllele, allele1, allele2, allele1Count, allele2Count, allele1AF, and allele2AF.

2. Run model likelihood calculation and selection.

RunModelLikelihood(
    seg = paste0(out_dir,"/",prefix,"_GATK_AI_segment.tsv"),
    out_dir = out_dir,
    prefix = prefix,
    gender = gender,
    modelminprobes = 20,
    modelminAIsize = 5000000,
    minsf = 0.4,
    callcov = 0.3,
    thread = 6)

Parameters for RunModelLikelihood

Notes:

• Parameters marked as (Optional) can be omitted and have defaults. - For a full description of all arguments and advanced options, see the function reference or ?RunModelLikelihood in R.

3. Run annotation segments.

AnnotateSegments(
    input = paste0(out_dir,"/",prefix,"_final_calls.tsv"),
    out_dir = out_dir,
    prefix = prefix,
    cytoband = cytoband,
    whitelist_edge = whitelist_edge,
    gene = gene)

Parameters for AnnotateSegments

4. Generating CNV plot

RunPlotCNV(
    seg = paste0(out_dir,"/",prefix,"_CNV_annotation.tsv"),
    cr =cr,
    ballele = ai,
    ai_binsize = 100000,
    cov_binsize = 100000,
    whitelist = whitelist_bed,
    gender = gender,
    out_dir = out_dir,
    prefix = prefix,
    aitype = "dragen"
  )

Parameters for RunPlotCNV

5. Generating AI segment quality file.

BafQC(
    annofile = paste0(out_dir,"/",prefix,"_CNV_annotation.tsv"),
    out_dir = out_dir,
    prefix = prefix)

Parameters for BafQC

Prepare input files

1. Run GATK in tumor-only mode by default parameters. Below is a summary of the GATK tumor-only mode command used in our pipeline. Please see the GATK website for details. Files will be used in XploR is sample.counts, sample.called.seg, sample.allelic_counts and sample.denoisedCR.tsv.
2. The allelic count file also could generate by other software like DRAGEN or samtools.

Supporting allelic count file format

Prepare Reference Files

Panel of normal reference

A Panel of Normals (PON) is required and should be generated using GATK, DRAGEN, or any other software capable of producing allelic count files.

Note: Male and female PON files need to be generated separately.

A. Whitelist, Blacklist, and Detectable Boundary Files

These files are generated from the PON HD5 file (from GATK), a cytoband file, and gender information. They are essential for downstream processing and include:

• Blacklist BED: Regions to exclude
• Whitelist BED: Regions to include
• Detectable Edge File: Defines detectable boundaries

These files are created based on the GATK Panel of Normals.
See the function documentation in R: ?PonProcess or help("PonProcess", package = "XploR").

Example usage:

PonProcess(
  pon_file = pon_hdh5_file,
  blacklist_bed = output_blacklist_bed,
  whitelist_bed = output_whitelist_bed,
  cytoband = cytoband,
  detectable_edge = output_detectable_edge,
  gender = gender
)

B.Panel of Normals Based on Allelic Count Files

The ai_pon_file should be a text file listing the paths to normal allelic count files generated by GATK, DRAGEN, or other software.

You can process these files to generate the PON reference for allelic imbalance using:

PONAIprocess(
  ai_pon_file = ai_pon_file,
  aitype = "GATK",
  minsnpcov = 20,
  output = "/Pathtoresults",
  prefix = "PONAI",
  maxgap = 2000000,
  maxbinsize = 5000000,
  minbinsize = 500000,
  snpnum = 30,
  gender = "female"
)

Parameters for PONAIprocess

Gene annotation reference:

Gene annotation can be obtained from various sources (e.g., Ensembl, UCSC, Gencode, RefSeq). An example file is included with the package:

gene <- system.file("extdata", "RefSeqCurated.genePred.gene_region.txt", package = "XploR")
head(read.table(gene, header = TRUE, sep = "\t"))

Cytoband annotation reference:

Cytoband annotation files are typically downloaded from UCSC. An example file is included:

cytoband <- system.file("extdata", "hg19_cytoBand.dat", package = "XploR")
head(read.table(cytoband, header = TRUE, sep = "\t"))

Altorithm

Binning Strategy for allelic count Data

The BinMaf function implements a flexible binning strategy for minor allele frequency (MAF) data, supporting both tumor samples and panels of normal (PoN) samples. The binning can be performed using either a fixed number of SNPs per bin with additional criteria to handle genomic gaps and bin size limits. Within each bin, Gaussian mixture modeling (GMM) is applied to identify clusters in the MAF distribution.

Key features:

• Flexible binning:
  Bins are created for each chromosome by grouping consecutive SNPs until one of the following conditions is met:
  • The number of SNPs in the bin reaches a specified target (e.g., 30 SNPs per bin).
  • The genomic span of the bin exceeds a maximum bin size (e.g., 2,000,000 bp).
  • The gap between consecutive SNPs exceeds a maximum allowed gap (e.g., 1,000,000 bp).
  • The bin size is at least the specified minimum size (e.g., 500,000 bp).

This strategy ensures that bins are of consistent size and SNP content, while avoiding the inclusion of widely separated SNPs in the same bin, and is robust for both tumor and normal samples.

Segmentation of MAF track

In addition to CBS (Circular Binary Segmentation), our pipeline supports a “merge” mode for segmentation based on minor allele frequency (MAF) values. While CBS is the default and recommended strategy, “merge” mode offers a step-wise, rule-based approach to combine adjacent MAF segments.

Step-wise Merging Strategy:

• The process begins with a minimal SNP counts, which is incrementally increased until reaching the final SNP count defined by “snpmin”.
• At each round, the following steps are performed:
  1. Remove Small Segments: Discard all segments smaller than the current minimum SNP count.
  2. Merge Adjacent Segments: Evaluate and merge adjacent segments if the difference in segment MAF values is less than or equal to the mergeai parameter (user-specified MAF difference threshold).

Note:

• CBS segmentation remains the default and is generally recommended for most workflows.
• The “merge” mode is useful for certain applications where a rule-based, stepwise merging of MAF segments is preferred.

MAF Bias Correction Using Panel of Normal (PoN) Allelic Counts

To correct for systemetic MAF bias estimates, we use a panel of normal (PoN) allelic count files as a reference. For each segment in the tumor or sample of interest, we compare the segment’s MAF to the distribution of MAF values observed in the PoN for the same genomic region. This process ensures that technical or locus-specific biases in MAF are removed only when the tumor segment does not show significant deviation from the normal reference.

Correction process: - For each segment, identify all overlapping PoN segments and extract their MAF values. - If the segment’s MAF is not significantly different from the PoN MAF distribution (assessed via a Wilcoxon test or a small absolute difference), apply a logit-based centering correction: - The segment MAF is transformed to the logit scale, centered by the PoN median MAF, and then inverse-logit transformed back and capped at 0.5. - If the segment’s MAF is significantly different from the PoN, the original segment MAF is retained (no correction is applied), thus preserving true biological signal. - This approach ensures that only technical or systematic biases are corrected, while real allelic imbalance events in the tumor are preserved.

This method uses the panel of normals as an adaptive reference, providing robust bias correction without shrinking true tumor signals.

Purity and diploid coverage scale factor estimation

Estimate a Beta–Binomial Over-Dispersion Parameter \(\theta\) from a Panel of Normals (PoN)

To accurately model over-dispersion in minor allele frequency (MAF) data, we estimate a beta-binomial dispersion parameter (\(\theta\)) using a panel of normal (PoN) samples. This allows us to account for extra-binomial variation and improves the likelihood calculation for each segment.

For each bin \(b\) and depth stratum:

• Observed variance across normals: \[ v_b = \mathrm{Var}_s(p_{sb}) \]
• Binomial expectation: \[ m_b = \mathbb{E}_s\left[\frac{p_{0b}(1-p_{0b})}{D_{sb}}\right] \approx p_{0b}(1-p_{0b}) \cdot \mathrm{mean}_s\left(\frac{1}{D_{sb}}\right) \]
• Representative depth in that (bin, stratum): \[ \tilde d_b = \mathrm{median}_s(D_{sb}) \]
• Moment estimator of \(\theta\) per (bin, stratum): \[ \widehat{\theta}_{b} = \max\left(\frac{v_b/m_b - 1}{\tilde d_b - 1},\ 0\right) \]

Within each depth stratum, we take a robust center (median) of \(\widehat{\theta}_b\) to obtain \(\theta\) for that stratum.

where:

• \(p_{sb}\): MAFs of per sample \(s\) per bin \(b\)
• \(p_{0b}\): Reference (PoN) mean MAF for bin \(b\)
• \(D_{sb}\): Median depth for sample \(s\) in bin \(b\)
• \(v_b\): Observed variance of MAF across samples in bin \(b\)
• \(m_b\): Expected binomial variance for bin \(b\)
• \(\tilde d_b\): Representative (median) depth in bin \(b\)
• \(\widehat{\theta}_b\): Estimated over-dispersion for bin \(b\)
• \(\theta\): Final over-dispersion parameter for the depth stratum

Prior assignment based on parsimony principle

Priors are assigned to each potential copy number combination based on the principle of parsimony, which favors simpler (biologically less complex) allele configurations. The biological difficulty level reflects the number of steps required to reach a given allele combination from the baseline diploid state (1,1), where each step represents either a gain or loss of one allele.

• Baseline State: The diploid state (1,1) is assigned a difficulty level of 1.
• Single-Step Changes: Combinations reachable from (1,1) in one step (single gain or loss) are assigned a difficulty level of 2.
• Two-Step Changes: Combinations requiring two steps (e.g., loss then gain, or two sequential gains) receive a difficulty level of 3.
• Three-, Four-, Five-Step Changes: Difficulty levels increase with the number of steps, up to 6 for the most complex.
• Special Considerations:
  • Whole Chromosome Duplication: Combinations like (2,2) are assigned higher difficulty due to their rarity.
  • Sequential Gains: Combinations with sequential gains of the same chromosome (e.g., (3,1)) are considered less difficult than those involving loss followed by gain (e.g., (2,0)). The prior for each allele combination is calculated using an exponential decay function controlled by a decay rate parameter λ. prior = exp(-λ × Bio_diff) where Bio_diff is the assigned biological difficulty score for the copy number configuration.

Tumor Copy Number Estimation

For each genomic segment, the model computes a range of possible tumor-specific copy numbers (CN_tumor) that could result from observed data under different cancer cell fractions (ccf):

• Tumor Copy Number Formula:

\[ CN_{tumor} = \frac{C_i \times 2 / (\mu \times 100) - (1 - \rho) \times 2 - \rho \times (1 - ccf) \times 2}{\rho \times ccf} \]

where:

• \(C_i\): Observed segment copy number

• \(\mu\): Diploid coverage scale factor

• \(\rho\): Tumor purity

• \(ccf\): Cancer cell fraction

• Combination Generation:
  For each potential tumor CN, all feasible major and minor allele combinations are generated and filtered based on biological plausibility.

• CCF Value Calculation:
  For non-diploid segments, CCF is calculated; for diploid segments, it is set as NA.

Likelihood Calculation

• MAF Likelihood:
  For each combination, the B allele frequency likelihood is computed using the beta distribution, parameterized by \(\alpha\) and \(\beta\), derived from the expected MAF:

\[ \mathrm{Beta}(\alpha, \beta) \]

where:

• \(\alpha = K \times \mathrm{BAF} + \epsilon\)

• \(\beta = K \times (1 - \mathrm{BAF}) + \epsilon\)

• \(K\): Beta-binomial precision parameter, estimated for each segment based on local read depth and the over-dispersion parameter (see below).

• \(\epsilon\): Small positive value for numerical stability

• Estimation of \(K\):

For each segment, \(K\) is calculated as:

\[ K = \frac{\text{depth}}{1 + (\text{depth} - 1) \cdot \theta} - 1 \]

where:

• \(\text{depth}\): Median SNP depth for the segment

• \(\theta\): Beta-binomial over-dispersion parameter, estimated from the panel of normals (PoN) for the corresponding depth stratum

• Posterior Likelihood:
  The posterior likelihood for each combination incorporates both the BAF likelihood and the prior, weighted by a factor \(\gamma\):

\[ \text{Posterior Likelihood} = \text{MAF Likelihood} \times (\text{Prior})^\gamma \]

where:

• \(\text{MAF Likelihood}\): Likelihood of the observed minor frequency under the current model
• \(\text{Prior}\): Prior probability assigned to the allele combination based on biological plausibility
• \(\gamma\): Weighting factor controlling the influence of the prior in the posterior calculation

Assigning Calls for Each Segment

The SelectCallpersegment() function refines and selects the most likely allele combinations for each genomic segment, handling both clonal and subclonal events, and incorporates coverage differences and prior knowledge.

• Initialization and Preprocessing:
  • Replace any zero BAF likelihoods with a minimum likelihood value to ensure all terms contribute meaningfully.
  • Compute expected coverage for each model and the difference (cov_diff) from observed coverage.
• Selection of Top Likelihood Models:
  • For each segment, identify the top two models with the highest MAF likelihoods.
  • If a subclonal event is likely (e.g., the second-ranked model has ccf > 0.3 and comparable likelihood), select it as a subclonal event.
  • If both major and minor copy numbers are equal, selection is based on cov_diff.
• Handling Models with minor = 0:
  • For segments with minor = 0, where MAF likelihood is unreliable, select the model with the smallest cov_diff.
  • Ensure consistency in selection for both minor = 0 and minor ≠ 0 cases.
• Post-Processing:
  • Calculate and store the log-transformed likelihood value (log_MAF_likelihood) for each selected model.

Output

Segmentation output

sample_GATK_AI_segment.tsv ( Generared by ?RunAIsegmentation function)

Raw likelihood results under each configuration

sample_likelihood_raw.tsv (Generated by ?RunModelLikelihood() function)

Allelic combiantion ressult with maximum likelihood under each configuration

sample_top_likelihood_calls.tsv ( Generated by ?SelectCallpersegment() function ) The format is simillar with sample_likelihood_raw.tsv, with best allelic combiantion is selected for each segment under each diploid coverage scale factor and tumor purity configuration.

Likelihood for each combination of diploid coverage scale factor and tumor purity

sample_Models_likelihood.tsv ( Generated by ?SelectFinalModel() function )

Final output of CNV calling

sample_final_calls.tsv (Generated by ?RunModelLikelihood() function)

Model selection plots

Likelihood dot plot: The plot displays the likelihood ranking for all combinations of diploid coverage scale factor and tumor purity. The vertical dashed line indicates the likelihood cutoff used to define Tier 1 models.

Model plot: The model plot displays the likelihood values of different models, which are calculated based on potential combinations of diploid coverage scale factor and tumor purity. In the plot, red indicates higher likelihood, while blue signifies lower likelihood. The light blue dot indicates the final model selected by XploR.

Tier1 Models Overall: This plot shows copy number calls for each combination of diploid coverage scale factor and tumor purity. Red indicates gain, blue indicates loss, and white indicates no change. Each configuration is labeled on the y-axis. By evaluating coverage and allelic imbalance patterns in this overview, you can identify the reasonable range of diploid coverage scale factors and tumor purity values. This helps guide reruns with optimized parameter ranges if needed.

Tier1 Models Zoom in: A zoomed-in view that makes the y-axis configurations more visible for detailed inspection.

QC Summary Table

sample_PASS_STAT_chr.txt ( Generated by ?BafQC() function )

Annotation file

sample_CNV_annotation.tsv ( Generated by ?AnnotateSegments() function, only unique columns are listed ).

ISCN calculation rules: 1. All segments will be reported with start and end cytoband in ISCN format. however certain considerations are made for the position of the centromere: a. In metacentric chromosomes, if a segment crosses the centromere and the gaps between the segment and the telomere on both sides are less than 5MB, only the chromosome number will be reported. b. In metacentric chromosomes, if a segment does not cross the centromere, and the gaps between the segment and the centromere and the telomere are both less than 5MB, the chromosome number followed by ‘p’ or ‘q’ will be reported. c. In acrocentric chromosomes, if the segment fulfills rule ‘b’ above, only the chromosome number will be reported.

CNV plot

sample_CNV_plot.png ( Generated by ?RunPlotCNV() function). The CNV Plot shows a genome-wide summary of the copy number (top track), B-allele frequency (BAF, second track) data, tumor fraction( ccf, third tract ) and quality of segment ( bottom track). The Copy Number (CN), on the Y-axis, is a linear count of the number of copies of each chromosome in the tumor cells, taking tumor purity and tumor fraction into account. Each chromosome is plotted as a set of dots that collectively show the estimated sequence coverage for the chromosome, and as a narrow turquoise line that shows the final CN call for the chromosome. The BAF plot shows the variant allele fraction of SNPs across the genome with the same coloration used in the Copy Number plot. When the copy number of a chromosome changes, the BAF plot for an affected chromosome splits due to imbalance in chromosome counts. The variance of B-allele frequencies is quite high so the splitting of the BAF may be difficult to discern. To assist with interpreting the BAF plot, a turquoise line is drawn at the median level to show the imbalance.

Model selection and rerun

Full function and parameter list

For a complete list of all functions and their parameters, please visit the XploR function reference.

Each function page includes detailed parameter descriptions, usage examples, and links to related documentation.